package com.ruoyi.common.utils.jkalman.jama;

import java.io.Serializable;

/**
 * @ClassName CholeskyDecomposition
 * @Description TODO
 * @Author jiayuan
 * @Date 2023/11/15 17:33
 * @Version 1.0
 */
public class CholeskyDecomposition implements Serializable {
    private double[][] L;
    private int n;
    private boolean isspd;

    public CholeskyDecomposition(Matrix Arg) {
        double[][] A = Arg.getArray();
        this.n = Arg.getRowDimension();
        this.L = new double[this.n][this.n];
        this.isspd = Arg.getColumnDimension() == this.n;

        for(int j = 0; j < this.n; ++j) {
            double[] Lrowj = this.L[j];
            double d = 0.0;

            int k;
            for(k = 0; k < j; ++k) {
                double[] Lrowk = this.L[k];
                double s = 0.0;

                for(int i = 0; i < k; ++i) {
                    s += Lrowk[i] * Lrowj[i];
                }

                Lrowj[k] = s = (A[j][k] - s) / this.L[k][k];
                d += s * s;
                this.isspd &= A[k][j] == A[j][k];
            }

            d = A[j][j] - d;
            this.isspd &= d > 0.0;
            this.L[j][j] = Math.sqrt(Math.max(d, 0.0));

            for(k = j + 1; k < this.n; ++k) {
                this.L[j][k] = 0.0;
            }
        }

    }

    public boolean isSPD() {
        return this.isspd;
    }

    public Matrix getL() {
        return new Matrix(this.L, this.n, this.n);
    }

    public Matrix solve(Matrix B) {
        if (B.getRowDimension() != this.n) {
            throw new IllegalArgumentException("Matrix row dimensions must agree.");
        } else if (!this.isspd) {
            throw new RuntimeException("Matrix is not symmetric positive definite.");
        } else {
            double[][] X = B.getArrayCopy();
            int nx = B.getColumnDimension();

            int k;
            int j;
            int i;
            for(k = 0; k < this.n; ++k) {
                for(j = 0; j < nx; ++j) {
                    for(i = 0; i < k; ++i) {
                        X[k][j] -= X[i][j] * this.L[k][i];
                    }

                    X[k][j] /= this.L[k][k];
                }
            }

            for(k = this.n - 1; k >= 0; --k) {
                for(j = 0; j < nx; ++j) {
                    for(i = k + 1; i < this.n; ++i) {
                        X[k][j] -= X[i][j] * this.L[i][k];
                    }

                    X[k][j] /= this.L[k][k];
                }
            }

            return new Matrix(X, this.n, nx);
        }
    }
}
